Create a new analysis directories.
- general directory
- stains directory
- for plots
- for output of summary results
- for baseline tables
* General packages...
* Genomic packages...
In this notebook we explore the files for the bulk RNAseq analyses and select the subset of data we need.
Here we obtain data from the CONVOCALS_downstream in plaques.
Genes.xlsx - list of genes of interest.Variants.xlsx - list of variant(s) of interest.library(openxlsx)
gene_list_df <- read.xlsx(paste0(TARGET_loc, "/targets.xlsx"), sheet = "Genes")
gene_list <- unlist(gene_list_df$Gene)
gene_list [1] "SCGB3A2" "LIX1" "IGSF9B" "ND6" "ALB" "OR10A4" "RASEF" "NEDD4" "TRIM49D1" "TCL1A" "ND4L" "ATP8" "FBXO15" "F5" "TMEM212"
[16] "PTPRD" "CYP46A1" "OR2T33" "SORBS2" "ITGA7" "RNASE1" "FOS" "HMOX1" "LAPTM5" "MMP9" "SIGLEC1" "FTL" "CD14" "HCST" "TIMP3"
[31] "CCL2" "SAT1" "CD163" "PTGDS" "LGALS9" "ACKR1" "NT5DC2" "TGFBI" "C1QC" "S100A9"
First we will load the data:
Now you need to choose which of the datasets you will use for the downstream analyses.
estimateSizeFactors()vst()Note: we advise using
vst()transformed normalized data
First, we extract the clinical data.
bulkRNA_meta_clin <- as.tibble(colData(AERNAScomboSE))
bulkRNA_meta_clin %<>%
# mutate(macrophages = factor(macrophages, levels = c("no staining", "minor staining", "moderate staining", "heavy staining"))) %>%
# mutate(smc = factor(smc, levels = c("no staining", "minor staining", "moderate staining", "heavy staining"))) %>%
# mutate(calcification = factor(calcification, levels = c("no staining", "minor staining", "moderate staining", "heavy staining"))) %>%
# mutate(collagen = factor(collagen, levels = c("no staining", "minor staining", "moderate staining", "heavy staining"))) %>%
# mutate(fat = factor(fat, levels = c("no fat", "< 40% fat", "> 40% fat"))) %>%
mutate(study_number_row = STUDY_NUMBER) %>%
as.data.frame() %>%
column_to_rownames("study_number_row")
names(bulkRNA_meta_clin)[names(bulkRNA_meta_clin) == "STUDY_NUMBER"] <- "study_number"
head(bulkRNA_meta_clin)[1] 1092 69
Second, we define the variables we need.
# Baseline table variables
basetable_vars = c("Hospital", "ORyear", "Artery_summary",
"Age", "Gender",
# "TC_finalCU", "LDL_finalCU", "HDL_finalCU", "TG_finalCU",
"TC_final", "LDL_final", "HDL_final", "TG_final",
# "hsCRP_plasma",
"systolic", "diastoli", "GFR_MDRD", "BMI",
"KDOQI", "BMI_WHO",
"SmokerStatus", "AlcoholUse",
"DiabetesStatus",
"Hypertension.selfreport", "Hypertension.selfreportdrug", "Hypertension.composite", "Hypertension.drugs",
"Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
"Stroke_Dx", "sympt", "Symptoms.5G", "AsymptSympt", "AsymptSympt2G",
"Symptoms.Update2G", "Symptoms.Update3G",
"restenos", "stenose", "Artery_summary",
"CAD_history", "PAOD", "Peripheral.interv",
"EP_composite", "EP_composite_time", "epcom.3years",
"EP_major", "EP_major_time","epmajor.3years",
"MAC_rankNorm", "SMC_rankNorm", "Macrophages.bin", "SMC.bin",
"Neutrophils_rankNorm", "MastCells_rankNorm",
"IPH.bin", "VesselDensity_rankNorm",
"Calc.bin", "Collagen.bin",
"Fat.bin_10", "Fat.bin_40",
"OverallPlaquePhenotype", "Plaque_Vulnerability_Index")
basetable_bin = c("Gender", "Artery_summary",
"KDOQI", "BMI_WHO",
"SmokerStatus", "AlcoholUse",
"DiabetesStatus",
"Hypertension.selfreport", "Hypertension.selfreportdrug", "Hypertension.composite", "Hypertension.drugs",
"Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
"Stroke_Dx", "sympt", "Symptoms.5G", "AsymptSympt", "AsymptSympt2G",
"Symptoms.Update2G", "Symptoms.Update3G",
"restenos", "stenose", "Artery_summary",
"CAD_history", "PAOD", "Peripheral.interv",
"EP_composite", "Macrophages.bin", "SMC.bin",
"IPH.bin",
"Calc.bin", "Collagen.bin",
"Fat.bin_10", "Fat.bin_40",
"OverallPlaquePhenotype", "Plaque_Vulnerability_Index")
# basetable_bin
basetable_con = basetable_vars[!basetable_vars %in% basetable_bin]
# basetable_conFrom here we can analyze whether specific genes differ between groups, or do this for the entire gene set as part of DE analysis, and then select our genes of interest.
From there, we extract the gene expression values, plus gene identifier, and select the genes of our interest.
expression_data <- assay(AERNAScomboSE)
# extract expression values from vsd, including ensembl names
expression_data <- as_tibble(data.frame(gene_ensembl = rowRanges(AERNAScomboSE)$feature_id, assay(AERNAScomboSE))) %>%
mutate_at(vars(c("gene_ensembl")), list(as.character)) ## gene_ensembl needs to be character for annotation to workmutate_at: no changes
# annotations
# gene symbol - via org.Hs.eg.db
# columns(org.Hs.eg.db)
expression_data$symbol <- mapIds(org.Hs.eg.db,
keys = expression_data$gene_ensembl,
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "first")'select()' returned 1:many mapping between keys and columns
# tidy and subset
expression_data_sel <- expression_data %>%
dplyr::select(gene_ensembl, symbol, everything()) %>%
# filter(symbol == "APOE" | symbol == "TRIB3") %>% # filter APOE and TRIB3
dplyr::filter(symbol %in% gene_list)
head(expression_data_sel)
# tidy and subset non-selected genes
set.seed(141619)
expression_data_sample <- expression_data %>%
dplyr::select(gene_ensembl, symbol, everything()) %>%
sample_n(1000) %>%
unite(symbol_ensembl, symbol, gene_ensembl, sep = "_", remove = FALSE)sample_n: removed 20,835 rows (95%), 1,000 rows remaining
expression_data_sample_mean <- expression_data_sample %>%
select_if(is.numeric) %>%
colMeans() %>%
as_tibble(rownames = "study_number") %>%
dplyr::rename(expression_value_sample = value)select_if: dropped 3 variables (symbol_ensembl, gene_ensembl, symbol)
Furthermore, the expression_data_sel df was gathered into a long form
df for annotation with symptoms variables from the
vsd object, and later visualization and statistics.
# gather expression_data_sel df into long df form for annotation, plotting and statistics
expression_long <-
gather(expression_data_sel, key = "study_number", value = "expression_value", -c(gene_ensembl, symbol))gather: reorganized (ae1, ae1026, ae1029, ae1032, ae1054, …) into (study_number, expression_value) [was 40x1094, now 43680x4]
# old school way
# Annotate with smoking variables
# sample_ids <- expression_long$study_number
# mm <- match(expression_long$study_number, rownames(colData(vsd)))
#
# ## Add traits to df
# ## Binary traits
# expression_long$sex <- colData(vsd)$sex[mm]
# expression_long$testosterone <- colData(vsd)$testosterone[mm]
# expression_long$t_e2_ratio <- colData(vsd)$t_e2_ratio[mm]
# new school way
plaque_phenotypes_cat <- c("Macrophages.bin",
"SMC.bin",
"Calc.bin",
"Collagen.bin",
"Fat.bin_10",
# "Fat.bin_40",
"IPH.bin")
plaque_phenotypes_num <- c("MAC_rankNorm", #"macmean0",
"SMC_rankNorm", #"smcmean0",
# "MastCells_rankNorm", #"mast_cells_plaque",
# "Neutrophils_rankNorm", #"neutrophils",
"VesselDensity_rankNorm")
expression_long <- expression_long %>%
left_join(bulkRNA_meta_clin %>% dplyr::select(study_number,
plaque_phenotypes_cat,
plaque_phenotypes_num,
epmajor.3years,
#epmajor.30days,
AsymptSympt2G,
Gender, Hospital, StudyName),
by = "study_number") %>%
mutate(epmajor_3years_yn = str_replace_all(epmajor.3years, c("Excluded" = "yes", "Included" = "no"))) #%>%left_join: added 14 columns (Macrophages.bin, SMC.bin, Calc.bin, Collagen.bin, Fat.bin_10, …) > rows only in x 0 > rows only in bulkRNA_meta_clin %>% d.. ( 0) > matched rows 43,680 > ======== > rows total 43,680
# mutate(epmajor.30days_yn = str_replace_all(epmajor.30days, c("Excluded" = "yes", "Included" = "no")))
head(expression_long)In case some genes are not available in our data we could filter them here.
[1] "SCGB3A2" "LIX1" "IGSF9B" "ND6" "ALB" "OR10A4" "RASEF" "NEDD4" "TRIM49D1" "TCL1A" "ND4L" "ATP8" "FBXO15" "F5" "TMEM212"
[16] "PTPRD" "CYP46A1" "OR2T33" "SORBS2" "ITGA7" "RNASE1" "FOS" "HMOX1" "LAPTM5" "MMP9" "SIGLEC1" "FTL" "CD14" "HCST" "TIMP3"
[31] "CCL2" "SAT1" "CD163" "PTGDS" "LGALS9" "ACKR1" "NT5DC2" "TGFBI" "C1QC" "S100A9"
This code is just an example to filter the list from genes that are not in the data.
Figure 1: Expression of genes of interest: boxplots
# Make directory for plots
ifelse(!dir.exists(file.path(QC_loc, "/Boxplots")),
dir.create(file.path(QC_loc, "/Boxplots")),
FALSE)[1] FALSE
BOX_loc = paste0(QC_loc,"/Boxplots")
for(GENE in gene_list_qc){
cat(paste0("Plotting expression for ", GENE,".\n"))
temp <- subset(expression_long, symbol == GENE)
compare_means(expression_value ~ Gender, data = temp)
p1 <- ggpubr::ggboxplot(temp,
x = "Gender",
y = "expression_value",
color = "Gender",
palette = "npg",
add = "jitter",
ylab = paste0("normalized expression ", GENE,"" ),
repel = TRUE
) + stat_compare_means()
print(p1)
cat(paste0("Saving image for ", GENE,".\n"))
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_gender.png"), plot = last_plot())
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_gender.pdf"), plot = last_plot())
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_gender.ps"), plot = last_plot())
compare_means(expression_value ~ Gender, data = temp)
p2 <- ggpubr::ggboxplot(temp,
x = "Gender",
y = "expression_value",
color = "StudyName",
palette = "npg",
add = "jitter",
ylab = paste0("normalized expression ", GENE,"" ),
repel = TRUE
) #+ stat_compare_means()
print(p2)
cat(paste0("Saving image for ", GENE,".\n"))
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_Gender_by_AERNAStudy.png"), plot = last_plot())
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_Gender_by_AERNAStudy.pdf"), plot = last_plot())
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_Gender_by_AERNAStudy.ps"), plot = last_plot())
compare_means(expression_value ~ StudyName, data = temp)
p3 <- ggpubr::ggboxplot(temp,
x = "StudyName",
y = "expression_value",
color = "StudyName",
palette = "npg",
add = "jitter",
ylab = paste0("normalized expression ", GENE,"" ),
repel = TRUE
) + stat_compare_means()
print(p3)
cat(paste0("Saving image for ", GENE,".\n"))
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy.png"), plot = last_plot())
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy.pdf"), plot = last_plot())
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy.ps"), plot = last_plot())
compare_means(expression_value ~ StudyName, data = temp)
p4 <- ggpubr::ggboxplot(temp,
x = "StudyName",
y = "expression_value",
color = "Gender",
palette = "npg",
add = "jitter",
ylab = paste0("normalized expression ", GENE,"" ),
repel = TRUE
) #+ stat_compare_means()
print(p4)
cat(paste0("Saving image for ", GENE,".\n"))
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy_by_Gender.png"), plot = last_plot())
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy_by_Gender.pdf"), plot = last_plot())
ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy_by_Gender.ps"), plot = last_plot())
rm(temp, p1, p2, p3, p4 )
}Plotting expression for SCGB3A2.
Saving image for SCGB3A2.
Saving image for SCGB3A2.
Saving image for SCGB3A2.
Saving image for SCGB3A2.
Plotting expression for LIX1.
Saving image for LIX1.
Saving image for LIX1.
Saving image for LIX1.
Saving image for LIX1.
Plotting expression for IGSF9B.
Saving image for IGSF9B.
Saving image for IGSF9B.
Saving image for IGSF9B.
Saving image for IGSF9B.
Plotting expression for ND6.
Saving image for ND6.
Saving image for ND6.
Saving image for ND6.
Saving image for ND6.
Plotting expression for ALB.
Saving image for ALB.
Saving image for ALB.
Saving image for ALB.
Saving image for ALB.
Plotting expression for OR10A4.
Saving image for OR10A4.
Saving image for OR10A4.
Saving image for OR10A4.
Saving image for OR10A4.
Plotting expression for RASEF.
Saving image for RASEF.
Saving image for RASEF.
Saving image for RASEF.
Saving image for RASEF.
Plotting expression for NEDD4.
Saving image for NEDD4.
Saving image for NEDD4.
Saving image for NEDD4.
Saving image for NEDD4.
Plotting expression for TRIM49D1.
Saving image for TRIM49D1.
Saving image for TRIM49D1.
Saving image for TRIM49D1.
Saving image for TRIM49D1.
Plotting expression for TCL1A.
Saving image for TCL1A.
Saving image for TCL1A.
Saving image for TCL1A.
Saving image for TCL1A.
Plotting expression for ND4L.
Saving image for ND4L.
Saving image for ND4L.
Saving image for ND4L.
Saving image for ND4L.
Plotting expression for ATP8.
Saving image for ATP8.
Saving image for ATP8.
Saving image for ATP8.
Saving image for ATP8.
Plotting expression for FBXO15.
Saving image for FBXO15.
Saving image for FBXO15.
Saving image for FBXO15.
Saving image for FBXO15.
Plotting expression for F5.
Saving image for F5.
Saving image for F5.
Saving image for F5.
Saving image for F5.
Plotting expression for TMEM212.
Saving image for TMEM212.
Saving image for TMEM212.
Saving image for TMEM212.
Saving image for TMEM212.
Plotting expression for PTPRD.
Saving image for PTPRD.
Saving image for PTPRD.
Saving image for PTPRD.
Saving image for PTPRD.
Plotting expression for CYP46A1.
Saving image for CYP46A1.
Saving image for CYP46A1.
Saving image for CYP46A1.
Saving image for CYP46A1.
Plotting expression for OR2T33.
Saving image for OR2T33.
Saving image for OR2T33.
Saving image for OR2T33.
Saving image for OR2T33.
Plotting expression for SORBS2.
Saving image for SORBS2.
Saving image for SORBS2.
Saving image for SORBS2.
Saving image for SORBS2.
Plotting expression for ITGA7.
Saving image for ITGA7.
Saving image for ITGA7.
Saving image for ITGA7.
Saving image for ITGA7.
Plotting expression for RNASE1.
Saving image for RNASE1.
Saving image for RNASE1.
Saving image for RNASE1.
Saving image for RNASE1.
Plotting expression for FOS.
Saving image for FOS.
Saving image for FOS.
Saving image for FOS.
Saving image for FOS.
Plotting expression for HMOX1.
Saving image for HMOX1.
Saving image for HMOX1.
Saving image for HMOX1.
Saving image for HMOX1.
Plotting expression for LAPTM5.
Saving image for LAPTM5.
Saving image for LAPTM5.
Saving image for LAPTM5.
Saving image for LAPTM5.
Plotting expression for MMP9.
Saving image for MMP9.
Saving image for MMP9.
Saving image for MMP9.
Saving image for MMP9.
Plotting expression for SIGLEC1.
Saving image for SIGLEC1.
Saving image for SIGLEC1.
Saving image for SIGLEC1.
Saving image for SIGLEC1.
Plotting expression for FTL.
Saving image for FTL.
Saving image for FTL.
Saving image for FTL.
Saving image for FTL.
Plotting expression for CD14.
Saving image for CD14.
Saving image for CD14.
Saving image for CD14.
Saving image for CD14.
Plotting expression for HCST.
Saving image for HCST.
Saving image for HCST.
Saving image for HCST.
Saving image for HCST.
Plotting expression for TIMP3.
Saving image for TIMP3.
Saving image for TIMP3.
Saving image for TIMP3.
Saving image for TIMP3.
Plotting expression for CCL2.
Saving image for CCL2.
Saving image for CCL2.
Saving image for CCL2.
Saving image for CCL2.
Plotting expression for SAT1.
Saving image for SAT1.
Saving image for SAT1.
Saving image for SAT1.
Saving image for SAT1.
Plotting expression for CD163.
Saving image for CD163.
Saving image for CD163.
Saving image for CD163.
Saving image for CD163.
Plotting expression for PTGDS.
Saving image for PTGDS.
Saving image for PTGDS.
Saving image for PTGDS.
Saving image for PTGDS.
Plotting expression for LGALS9.
Saving image for LGALS9.
Saving image for LGALS9.
Saving image for LGALS9.
Saving image for LGALS9.
Plotting expression for ACKR1.
Saving image for ACKR1.
Saving image for ACKR1.
Saving image for ACKR1.
Saving image for ACKR1.
Plotting expression for NT5DC2.
Saving image for NT5DC2.
Saving image for NT5DC2.
Saving image for NT5DC2.
Saving image for NT5DC2.
Plotting expression for TGFBI.
Saving image for TGFBI.
Saving image for TGFBI.
Saving image for TGFBI.
Saving image for TGFBI.
Plotting expression for C1QC.
Saving image for C1QC.
Saving image for C1QC.
Saving image for C1QC.
Saving image for C1QC.
Plotting expression for S100A9.
Saving image for S100A9.
Saving image for S100A9.
Saving image for S100A9.
Saving image for S100A9.
Figure 2A: Expression of genes of interest: histograms
# Make directory for plots
ifelse(!dir.exists(file.path(QC_loc, "/Histograms")),
dir.create(file.path(QC_loc, "/Histograms")),
FALSE)[1] FALSE
HISTOGRAM_loc = paste0(QC_loc,"/Histograms")
for(GENE in gene_list_qc){
# cat(paste0("Plotting expression for ", GENE,".\n"))
temp <- subset(expression_long, symbol == GENE)
p1 <- ggpubr::gghistogram(temp,
x = "expression_value",
y = "..count..",
color = "Gender", fill = "Gender",
palette = "npg",
add = "median",
ylab = paste0("normalized expression ", GENE,"" )
)
print(p1)
cat(paste0("Saving image for ", GENE,".\n"))
ggsave(filename = paste0(HISTOGRAM_loc, "/", Today, ".",GENE,".distribution.png"), plot = last_plot())
ggsave(filename = paste0(HISTOGRAM_loc, "/", Today, ".",GENE,".distribution.pdf"), plot = last_plot())
ggsave(filename = paste0(HISTOGRAM_loc, "/", Today, ".",GENE,".distribution.ps"), plot = last_plot())
rm(temp, p1 )
}Saving image for SCGB3A2.
Saving image for LIX1.
Saving image for IGSF9B.
Saving image for ND6.
Saving image for ALB.
Saving image for OR10A4.
Saving image for RASEF.
Saving image for NEDD4.
Saving image for TRIM49D1.
Saving image for TCL1A.
Saving image for ND4L.
Saving image for ATP8.
Saving image for FBXO15.
Saving image for F5.
Saving image for TMEM212.
Saving image for PTPRD.
Saving image for CYP46A1.
Saving image for OR2T33.
Saving image for SORBS2.
Saving image for ITGA7.
Saving image for RNASE1.
Saving image for FOS.
Saving image for HMOX1.
Saving image for LAPTM5.
Saving image for MMP9.
Saving image for SIGLEC1.
Saving image for FTL.
Saving image for CD14.
Saving image for HCST.
Saving image for TIMP3.
Saving image for CCL2.
Saving image for SAT1.
Saving image for CD163.
Saving image for PTGDS.
Saving image for LGALS9.
Saving image for ACKR1.
Saving image for NT5DC2.
Saving image for TGFBI.
Saving image for C1QC.
Saving image for S100A9.
Figure 2B: Expression of genes of interest: density plots
# Make directory for plots
ifelse(!dir.exists(file.path(QC_loc, "/Density")),
dir.create(file.path(QC_loc, "/Density")),
FALSE)[1] FALSE
DENSITY_loc = paste0(QC_loc,"/Density")
for(GENE in gene_list_qc){
# cat(paste0("Plotting expression for ", GENE,".\n"))
temp <- subset(expression_long, symbol == GENE)
p1 <- ggpubr::gghistogram(temp,
x = "expression_value",
y = "..density..",
color = "Gender", fill = "Gender",
palette = "npg",
add = "median",
ylab = paste0("normalized expression ", GENE,"" )
)
print(p1)
cat(paste0("Saving image for ", GENE,".\n"))
ggsave(filename = paste0(DENSITY_loc, "/", Today, ".",GENE,".density.png"), plot = last_plot())
ggsave(filename = paste0(DENSITY_loc, "/", Today, ".",GENE,".density.pdf"), plot = last_plot())
ggsave(filename = paste0(DENSITY_loc, "/", Today, ".",GENE,".density.ps"), plot = last_plot())
rm(temp, p1 )
}Saving image for SCGB3A2.
Saving image for LIX1.
Saving image for IGSF9B.
Saving image for ND6.
Saving image for ALB.
Saving image for OR10A4.
Saving image for RASEF.
Saving image for NEDD4.
Saving image for TRIM49D1.
Saving image for TCL1A.
Saving image for ND4L.
Saving image for ATP8.
Saving image for FBXO15.
Saving image for F5.
Saving image for TMEM212.
Saving image for PTPRD.
Saving image for CYP46A1.
Saving image for OR2T33.
Saving image for SORBS2.
Saving image for ITGA7.
Saving image for RNASE1.
Saving image for FOS.
Saving image for HMOX1.
Saving image for LAPTM5.
Saving image for MMP9.
Saving image for SIGLEC1.
Saving image for FTL.
Saving image for CD14.
Saving image for HCST.
Saving image for TIMP3.
Saving image for CCL2.
Saving image for SAT1.
Saving image for CD163.
Saving image for PTGDS.
Saving image for LGALS9.
Saving image for ACKR1.
Saving image for NT5DC2.
Saving image for TGFBI.
Saving image for C1QC.
Saving image for S100A9.
Figure 3: comparing expression of genes of interest to mean expression of a sample of 1,000 random genes
expression_wide <- expression_long %>%
tidyr::unite("symbol_ensembl", symbol:gene_ensembl, remove = TRUE, sep = "_") %>%
# dplyr::select(-gene_ensembl, -symbol) %>%
tidyr::spread(key = symbol_ensembl, value = expression_value) # the next 3 lines of code gave an error when selecting for genes_interest, since one of the genes of interest is missing: FGF3 is not in the data set. So, we need to select for the other 15 genes.
temp <- expression_long %>%
tidyr::unite("symbol_ensembl", symbol:gene_ensembl, remove = TRUE, sep = "_") %>%
dplyr::select(symbol_ensembl)
temp <- temp %>%
tidyr::separate(symbol_ensembl, c("symbol", "gene_ensembl"), remove = FALSE) %>%
dplyr::distinct(., symbol_ensembl, .keep_all = TRUE)
symbol_ensembl_gene_list_qc <- temp[temp$symbol %in% gene_list_qc,]$symbol_ensemblexpression_wide2 <- expression_wide %>%
left_join(expression_data_sample_mean, by = "study_number") %>%
dplyr::select(study_number, StudyName, Gender, symbol_ensembl_gene_list_qc, expression_value_sample)left_join: added one column (expression_value_sample) > rows only in x 0 > rows only in expression_data_sample_.. ( 0) > matched rows 1,092 > ======= > rows total 1,092
expression_long2 <- expression_wide2 %>%
gather(gene, expression_value, -study_number, -StudyName, -Gender) %>%
mutate(gene = str_replace_all(gene, c("expression_value_sample" = "Random genes"))) %>%
tidyr::separate(gene, c("gene", "symbol_ensembl"), remove = FALSE, sep = "_") %>%
dplyr::select(-symbol_ensembl) %>%
dplyr::mutate(gene = factor(gene, levels = c("Random genes", gene_list_qc)))gather: reorganized (CYP46A1_ENSG00000036530, NEDD4_ENSG00000069869, IGSF9B_ENSG00000080854, FTL_ENSG00000087086, SIGLEC1_ENSG00000088827, …) into (gene, expression_value) [was 1092x44, now 44772x5]Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1092 rows [43681, 43682, 43683, 43684, 43685, 43686, 43687, 43688, 43689, 43690, 43691, 43692, 43693, 43694, 43695, 43696, 43697, 43698, 43699, 43700, ...].
mean_1000_genes <- mean(expression_data_sample_mean$expression_value_sample)
# head(expression_long2)
#p1 <- ggpubr::ggboxplot(expression_long2,
x = "gene",
y = "expression_value",
color = uithof_color[16],
add = "jitter",
add.params = list(size = 3, jitter = 0.3),
ylab = paste0("normalized expression ")
) +
geom_hline(yintercept = mean_1000_genes, linetype = "dashed", color = uithof_color[26], size = 0.4) +
theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), # change orientation of x-axis labels
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14))
p1
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes.png"), plot = last_plot())Saving 7.29 x 4.51 in image
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes.pdf"), plot = last_plot())Saving 7.29 x 4.51 in image
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes.ps"), plot = last_plot())Saving 7.29 x 4.51 in image
p2 <- ggpubr::ggboxplot(expression_long2,
x = "gene",
y = "expression_value",
color = "StudyName",
palette = "npg",
add = "jitter",
add.params = list(size = 3, jitter = 0.3),
ylab = paste0("normalized expression ")
) +
geom_hline(yintercept = mean_1000_genes, linetype = "dashed", color = uithof_color[26], size = 0.4) +
theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), # change orientation of x-axis labels
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14))
p2
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_AERNAStudy.png"), plot = last_plot())Saving 7.29 x 4.51 in image
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_AERNAStudy.pdf"), plot = last_plot())Saving 7.29 x 4.51 in image
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_AERNAStudy.ps"), plot = last_plot())Saving 7.29 x 4.51 in image
p3 <- ggpubr::ggboxplot(expression_long2,
x = "gene",
y = "expression_value",
color = "Gender",
palette = "npg",
add = "jitter",
add.params = list(size = 3, jitter = 0.3),
ylab = paste0("normalized expression ")
) +
geom_hline(yintercept = mean_1000_genes, linetype = "dashed", color = uithof_color[26], size = 0.4) +
theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), # change orientation of x-axis labels
axis.text.y = element_text(size = 12),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14))
p3
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_Gender.png"), plot = last_plot())Saving 7.29 x 4.51 in image
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_Gender.pdf"), plot = last_plot())Saving 7.29 x 4.51 in image
ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_Gender.ps"), plot = last_plot())Saving 7.29 x 4.51 in image
List all the gene for all unique genes where their mean
expression is higher than the mean of the 1,000 random genes
(mean_1000_genes).
# calculate the mean for each gene in expression_long2
expression_long2_mean <- expression_long2 %>%
group_by(gene) %>%
summarise(mean_expression = mean(expression_value))summarise: now 41 rows and 2 columns, ungrouped
expression_long2_high <- expression_long2_mean %>%
dplyr::filter(mean_expression > mean_1000_genes) %>%
dplyr::select(gene) %>%
dplyr::distinct(gene, .keep_all = TRUE)
expression_long2_highList all genes with expression lower than
mean_1000_genes.
expression_long2_low <- expression_long2_mean %>%
dplyr::filter(mean_expression < mean_1000_genes) %>%
dplyr::select(gene) %>%
dplyr::distinct(gene, .keep_all = TRUE)
expression_long2_lowList all genes where the mean of their expression is higher or lower than the mean of 1,000 random genes, and where there is a statistical difference between gender.
# 1. Calculate mean expression per gender per gene
expression_long2_mean_gender <- expression_long2 %>%
dplyr::group_by(gene, Gender) %>%
dplyr::summarise(mean_expression = mean(expression_value), .groups = 'drop')
# Print the mean expression per gender
print(expression_long2_mean_gender)
# 2. Perform a Wilcoxon rank-sum test (Mann-Whitney U test) for gender differences per gene
expression_long2_mean_gender_t_test_results <- expression_long2 %>%
dplyr::group_by(gene) %>%
dplyr::summarise(wilcox_test = list(wilcox.test(expression_value ~ Gender)), .groups = 'drop') %>%
dplyr::mutate(p_value = sapply(wilcox_test, function(x) x$p.value)) %>%
dplyr::filter(p_value <= 0.05) %>% # Filter for p-values below 0.05
dplyr::arrange(p_value) # Order by p-value in ascending order
# Print the p-values from the t-tests
print(expression_long2_mean_gender_t_test_results)
# 3. Write the results to an Excel file using openxlsx
# Create a new workbook
wb <- createWorkbook()
# Add a worksheet
addWorksheet(wb, "Significant Genes")
# Write the dataframe to the worksheet
writeData(wb, "Significant Genes", expression_long2_mean_gender_t_test_results)
# Save the workbook to an Excel file
saveWorkbook(wb, paste0(OUT_loc, "/", Today,".iph_area_075_vst_dgea_gender_significant_genes.xlsx"), overwrite = TRUE)If we would put these correlations in one simple and comprehensible figure, we could use a correlation heatmap. Again, correlation coefficients used here are Spearman’s.
Figure 4: correlation heatmap between expression of genes of interest
library(tidyverse)
library(magrittr)
temp <- expression_wide %>%
column_to_rownames("study_number") %>%
dplyr::select(symbol_ensembl_gene_list_qc) %>%
drop_na() %>% # drop NA
Filter(function(x) sd(x) != 0, .) # filter variables with sd = 0
# average out columns that start with the same gene symbol, where the gene symbol is the first part of the column name separated by '_'; rename the columns to the part before '_'
# Load necessary library
library(dplyr)
# Assume your data frame is named `expression_wide`
# Step 1: Identify the columns that correspond to gene symbols using regex
gene_columns <- grep("ENSG\\d+", colnames(temp), value = TRUE)
# Step 2: Extract gene symbols using regular expression
# Extract everything before the underscore (_) as the gene symbol
gene_symbols <- sub("(_ENSG\\d+)", "", gene_columns)
# Step 3: Create a list of unique gene symbols
unique_genes <- unique(gene_symbols)
# Step 4: For each unique gene, either average the columns (if multiple) or rename (if single)
temp2 <- temp
for (gene in unique_genes) {
# Find all columns that match the gene symbol
gene_cols <- grep(paste0("^", gene, "_ENSG\\d+"), colnames(temp2), value = TRUE)
if (length(gene_cols) > 1) {
# If multiple columns for the gene, compute the row-wise average
temp2[[gene]] <- rowMeans(temp2[gene_cols], na.rm = TRUE)
# Remove the original ENSG columns after averaging
temp2 <- temp2 %>% dplyr::select(-all_of(gene_cols))
} else if (length(gene_cols) == 1) {
# If only one column for the gene, just rename it by removing the ENSG part
colnames(temp2)[colnames(temp2) == gene_cols] <- gene
}
}
# View the result
print(temp2)p1 <- pheatmap(data.matrix(temp.cor),
scale = "none",
cluster_rows = TRUE,
cluster_cols = TRUE,
legend = TRUE,
fontsize = 10)
p1We are saving the final list of genes of interest
We will create a list of samples that should be included based on
CEA, and having the proper informed consent (‘academic’). We will save
the SummarizedExperiment as a RDS file for easy loading.
The clinical data will also be saved as a separate
txt-file.
To this end you need to lookup the EnsemblID (see above) for your target(s).
# targets = c("ENSG00000166033", # target
# "ENSG00000169174") # positive control PCSK9
# targets = c("ENSG00000169174") # target/positive control PCSK9
targets = unlist(temp$gene_ensembl)
rm(temp)We make a SummarizedExperiment for the subset of RNAseq
data. We make sure to only include the samples we need based on informed
consent and we include only the requested variables.
Next, we are constructing the SummarizedExperiment.
AERNAS_meta_data = c("SampleType", "RNAseqTech", "RNAseqType", "RNAseqQC",
"StudyType", "StudyName", "Gender", "StudyBiobank", "SampleSize")
cat("\n* making a SummarizedExperiment ...\n")
* making a SummarizedExperiment ...
# this is all the data passing RNAseq quality control and UMI-corrected
# - after filtering on informed consent and artery type, the end sample size should be 1093
# - after filtering on 'no commercial business' based on informed consent, there are fewer samples: yyyy
cat(" > getting counts\n") > getting counts
[1] 40 1092
> meta data and clinical data
# bulkRNA_meta_clin_COMMERCIAL <- subset(bulkRNA_meta_clin, select = c("study_number", basetable_vars))
bulkRNA_meta_clin_ACADEMIC <- subset(bulkRNA_meta_clin, select = c("study_number",
AERNAS_meta_data, # meta data of RNAseq experiment
basetable_vars)) # selection of variables
dim(bulkRNA_meta_clin_ACADEMIC)[1] 1092 68
> rowRanges
# bulkRNA_meta_clin_COMMERCIAL <- subset(bulkRNA_meta_clin, select = c("study_number", basetable_vars))
bulkRNA_rowRanges <- rowRanges(AERNAScomboSE[targets])
bulkRNA_rowRangesGRanges object with 40 ranges and 2 metadata columns:
seqnames ranges strand | feature_id symbol
<Rle> <IRanges> <Rle> | <character> <character>
ENSG00000036530 14 99684298-99727301 + | ENSG00000036530 CYP46A1
ENSG00000069869 15 55826922-55993660 - | ENSG00000069869 NEDD4
ENSG00000080854 11 133896438-133956968 - | ENSG00000080854 IGSF9B
ENSG00000087086 19 48965309-48966879 + | ENSG00000087086 FTL
ENSG00000088827 20 3686970-3712600 - | ENSG00000088827 SIGLEC1
... ... ... ... . ... ...
ENSG00000198734 1 169511951-169586588 - | ENSG00000198734 F5
ENSG00000212907 MT 10470-10766 + | ENSG00000212907 ND4L
ENSG00000213088 1 159203307-159206500 + | ENSG00000213088 ACKR1
ENSG00000223417 11 89911111-89922245 - | ENSG00000223417 TRIM49D1
ENSG00000228253 MT 8366-8572 + | ENSG00000228253 ATP8
-------
seqinfo: 331 sequences from an unspecified genome; no seqlengths
> construction of the SE
(AERNAScomboSE_target <- SummarizedExperiment(assays = list(counts = as.matrix(bulkRNA_countsFilt)),
colData = bulkRNA_meta_clin_ACADEMIC,
rowRanges = bulkRNA_rowRanges,
metadata = "Athero-Express Biobank Study bulk RNA sequencing. Sample type: carotid plaques. Technology: CEL2-seq adapted for bulk RNA sequencing, thus 3'-focused. UMI-corrected. VST and size factor normalized. Target subset."))class: RangedSummarizedExperiment
dim: 40 1092
metadata(1): ''
assays(1): counts
rownames(40): ENSG00000036530 ENSG00000069869 ... ENSG00000223417 ENSG00000228253
rowData names(2): feature_id symbol
colnames(1092): ae1 ae1026 ... ae986 ae992
colData names(68): study_number SampleType ... OverallPlaquePhenotype Plaque_Vulnerability_Index
* removing intermediate files ...
Sample list and clinical data.
temp <- as.tibble(subset(colData(AERNAScomboSE_target), select = c("study_number", AERNAS_meta_data)))
fwrite(temp,
file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.samplelist.after_qc.IC_academic.csv"),
sep = ",", row.names = FALSE, col.names = TRUE,
showProgress = TRUE)
rm(temp)
temp <- as.tibble(subset(colData(AERNAScomboSE_target), select = c("study_number", basetable_vars)))
fwrite(temp,
file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.clinicaldata.after_qc.IC_academic.csv"),
sep = ",", row.names = FALSE, col.names = TRUE,
showProgress = TRUE)
rm(temp)Assay data and gene information.
temp <- as_tibble(assay(AERNAScomboSE_target))
fwrite(temp,
file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.assay.after_qc.IC_academic.csv"),
sep = ",", row.names = FALSE, col.names = TRUE,
showProgress = TRUE)
rm(temp)
temp <- as_tibble(rowRanges(AERNAScomboSE_target))
fwrite(temp,
file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.genedata.after_qc.IC_academic.csv"),
sep = ",", row.names = FALSE, col.names = TRUE,
showProgress = TRUE)
rm(temp)The subset as SummarizedExperiment() dataset in
RDS-format.
Here we just do a sanity check and compare the expression for a favorite gene.
ggpubr::gghistogram(as.tibble(t(subset(assay(AERNAScomboSE_target), AERNAScomboSE_target@rowRanges$symbol == "HMOX1"))),
x = "ENSG00000100292",
xlab = "HMOX1 (ENSG00000100292) expression\nlog-transformed, size corrected counts",
color = "white", fill = uithof_color[20],
rug = F, add_density = F,
add = c("median"),
add.params = list(color = uithof_color[3]),
ggtheme = theme_pubclean())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.Warning: `geom_vline()`: Ignoring `data` because `xintercept` was provided.
ggpubr::gghistogram(as.tibble(t(subset(assay(AERNAScomboSE_target), AERNAScomboSE_target@rowRanges$symbol == "SORBS2"))),
x = "ENSG00000154556",
xlab = "SORBS2 (ENSG00000154556) expression\nlog-transformed, size corrected counts",
color = "white", fill = uithof_color[16],
rug = F, add_density = F,
add = c("median"),
add.params = list(color = uithof_color[3]),
ggtheme = theme_pubclean())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.Warning: `geom_vline()`: Ignoring `data` because `xintercept` was provided.
Version: v1.1.0
Last update: 2024-10-17
Written by: Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
Description: Script to load bulk RNA sequencing data, and perform gene expression analyses, and visualisations.
Minimum requirements: R version 3.5.2 (2018-12-20) -- 'Eggshell Igloo', macOS Mojave (10.14.2).
**MoSCoW To-Do List**
The things we Must, Should, Could, and Would have given the time we have.
_M_
_S_
_C_
_W_
**Changes log**
* v1.1.0 Added more detailed sub-analyses on the mean expression by gender.
* v1.0.0 Inital version. Re-arranged to explore the subset of data.
R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS 15.1
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats4 grid tools stats graphics grDevices utils datasets methods base
other attached packages:
[1] magrittr_2.0.3 rmarkdown_2.28 annotables_0.2.0 EnhancedVolcano_1.22.0
[5] ggrepel_0.9.6 EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.28.1 AnnotationFilter_1.28.0
[9] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 mygene_1.40.0 txdbmaker_1.0.1 org.Hs.eg.db_3.19.1
[13] DESeq2_1.44.0 SummarizedExperiment_1.34.0 MatrixGenerics_1.16.0 matrixStats_1.4.1
[17] GenomicFeatures_1.56.0 GenomicRanges_1.56.2 AnnotationDbi_1.66.0 Biobase_2.64.0
[21] GenomeInfoDb_1.40.1 IRanges_2.38.1 S4Vectors_0.42.1 BiocGenerics_0.50.0
[25] Hmisc_5.1-3 survminer_0.4.9 survival_3.7-0 GGally_2.2.1
[29] PerformanceAnalytics_2.0.4 xts_0.14.0 zoo_1.8-12 ggcorrplot_0.1.4.999
[33] corrr_0.4.4 reshape2_1.4.4 meta_7.0-0 metadat_1.2-0
[37] qqman_0.1.9 sjlabelled_1.2.0 patchwork_1.3.0.9000 labelled_2.13.0
[41] sjPlot_2.8.16 UpSetR_1.4.0 ggpubr_0.6.0.999 forestplot_3.1.5
[45] abind_1.4-8 checkmate_2.3.2 pheatmap_1.0.12 devtools_2.4.5
[49] usethis_3.0.0 BlandAltmanLeh_0.3.1 tableone_0.13.2 openxlsx_4.2.7.1
[53] haven_2.5.4 eeptools_1.2.5 DT_0.33 knitr_1.48
[57] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 purrr_1.0.2
[61] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0 data.table_1.16.2
[65] naniar_1.1.0 tidylog_1.1.0 tidyr_1.3.1 dplyr_1.1.4
[69] optparse_1.7.5 readr_2.1.5 pander_0.6.5 R.utils_2.12.3
[73] R.oo_1.26.0 R.methodsS3_1.8.2 worcs_0.1.15 credentials_2.0.2
loaded via a namespace (and not attached):
[1] sjmisc_2.8.10 progress_1.2.3 urlchecker_1.0.1 nnet_7.3-19 Biostrings_2.72.1 vctrs_0.6.5
[7] rticles_0.27 digest_0.6.37 png_0.1-8 proxy_0.4-27 renv_1.0.11 MASS_7.3-61
[13] httpuv_1.6.15 withr_3.0.1 xfun_0.48 ellipsis_0.3.2 memoise_2.0.1 profvis_0.4.0
[19] ggsci_3.2.0 systemfonts_1.1.0 calibrate_1.7.7 ragg_1.3.3 Formula_1.2-5 sys_3.4.3
[25] prettyunits_1.2.0 datawizard_0.13.0 KEGGREST_1.44.1 promises_1.3.0 httr_1.4.7 rstatix_0.7.2
[31] restfulr_0.0.15 rstudioapi_0.16.0 UCSC.utils_1.0.0 miniUI_0.1.1.1 generics_0.1.3 base64enc_0.1-3
[37] curl_5.2.3 mitools_2.4 zlibbioc_1.50.0 ggeffects_1.7.2 GenomeInfoDbData_1.2.12 quadprog_1.5-8
[43] SparseArray_1.4.8 xtable_1.8-4 evaluate_1.0.1 S4Arrays_1.4.1 BiocFileCache_2.12.0 hms_1.1.3
[49] colorspace_2.1-1 filelock_1.0.3 getopt_1.20.4 lmtest_0.9-40 later_1.3.2 lattice_0.22-6
[55] XML_3.99-0.17 survey_4.4-2 class_7.3-22 pillar_1.9.0 nlme_3.1-166 performance_0.12.3
[61] compiler_4.4.1 stringi_1.8.4 minqa_1.2.8 GenomicAlignments_1.40.0 plyr_1.8.9 crayon_1.5.3
[67] BiocIO_1.14.0 chron_2.3-61 locfit_1.5-9.10 bit_4.5.0 mathjaxr_1.6-0 textshaping_0.4.0
[73] codetools_0.2-20 clisymbols_1.2.0 openssl_2.2.2 crosstalk_1.2.1 bslib_0.8.0 e1071_1.7-16
[79] mime_0.12 splines_4.4.1 Rcpp_1.0.13 dbplyr_2.5.0 blob_1.2.4 utf8_1.2.4
[85] lme4_1.1-35.5 fs_1.6.4 prereg_0.6.0 pkgbuild_1.4.4 gh_1.4.1 ggsignif_0.6.4
[91] sqldf_0.4-11 Matrix_1.7-0 tzdb_0.4.0 visdat_0.6.0 pkgconfig_2.0.3 cachem_1.1.0
[97] RSQLite_2.3.7 DBI_1.2.3 numDeriv_2016.8-1.1 fastmap_1.2.0 scales_1.3.0 Rsamtools_2.20.0
[103] broom_1.0.7 sass_0.4.9 coda_0.19-4.1 ggstats_0.7.0 insight_0.20.5 carData_3.0-5
[109] rpart_4.1.23 farver_2.1.2 gsubfn_0.7 yaml_2.3.10 foreign_0.8-87 sjstats_0.19.0
[115] rtracklayer_1.64.0 cli_3.6.3 lifecycle_1.0.4 askpass_1.2.1 sessioninfo_1.2.2 backports_1.5.0
[121] BiocParallel_1.38.0 timechange_0.3.0 gtable_0.3.5 rjson_0.2.23 metafor_4.6-0 parallel_4.4.1
[127] jsonlite_1.8.9 bitops_1.0-9 bit64_4.5.2 proto_1.0.0 zip_2.3.1 ranger_0.16.0
[133] jquerylib_0.1.4 survMisc_0.5.6 lazyeval_0.2.2 shiny_1.9.1 htmltools_0.5.8.1 KMsurv_0.1-5
[139] rappdirs_0.3.3 tinytex_0.53 glue_1.8.0 httr2_1.0.5 XVector_0.44.0 RCurl_1.98-1.16
[145] gridExtra_2.3 boot_1.3-31 R6_2.5.1 arm_1.14-4 km.ci_0.5-6 vcd_1.4-13
[151] CompQuadForm_1.4.3 labeling_0.4.3 cluster_2.1.6 pkgload_1.4.0 nloptr_2.1.1 DelayedArray_0.30.1
[157] tidyselect_1.2.1 ProtGenerics_1.36.0 htmlTable_2.4.3 xml2_1.3.6 car_3.1-3 munsell_0.5.1
[163] htmlwidgets_1.6.4 RColorBrewer_1.1-3 biomaRt_2.60.1 rlang_1.1.4 gert_2.1.4 remotes_2.5.0
[169] fansi_1.0.6
| © 1979-2024 Sander W. van der Laan | s.w.vanderlaan[at]gmail[dot]com | vanderlaanand.science. |